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ABSTRACT 


We present the results obtained from linear stability analysis and 2.5- 
dimensional magnetohydro dynamic (MHD) simulations of the magnetorotational 
instability (MRI), including the effects of cosmic rays (CRs). We took into ac¬ 
count of the CR diffusion along the magnetic held but neglect the cross-held-line 
diffusion. Two models are considered in this paper: shearing box model and 
differentially rotating cylinder model. We studied how MRI is affected by the 
initial CR pressure (i.e., energy) distribution. In the shearing box model, the 
initial state is uniform distribution. Linear analysis shows that the growth rate 
of MRI does not depend on the value of CR diffusion coefficient. In the differen¬ 
tially rotating cylinder model, the initial state is a constant angular momentum 
polytropic disk threaded by weak uniform vertical magnetic held. Linear analysis 
shows that the growth rate of MRI becomes larger if the CR diffusion coefficient 
is larger. Both results are confirmed by MHD simulations. The MHD simulation 
results show that the outward movement of matter by the growth of MRI is not 
impeded by the CR pressure gradient, and the centrifugal force which acts to the 
concentrated matter becomes larger. Consequently, the growth rate of MRI is in¬ 
creased. On the other hand, if the initial CR pressure is uniform, then the growth 
rate of the MRI barely depends on the value of the CR diffusion coefficient. 

Subject headings: accretion, accretion disks - cosmic rays - diffusion - instabilities - 
magnetic fields - MHD - Galaxy: disk 



3 


INTRODUCTION 


Magnetic field, an important component of the interstellar medium (ISM), is thought to 
be a key player in various active astrophysical phenomena. However, the dynamical role of 
cosmic ray (CR) (another component of ISM) in astrophysical activities has been underrated 


for quite a long time, although the energy density_of C 


magnetic field and turbulent gas motions (e.g., 


Parker 


Ls is of the same order as 


1969 


Ferriere 


Ginzburg & Ptuskin 


hat of 


1976 


20011 ). Still some effort have been made over the years. The most convenient way 


to study the effect of CRs on plasma flow is to desc ribe the system as a multi-fluid system 


where plasma and CRs are considered as fluid s (e.g., 


WebbetaL 


1986 


Webb 


1987|; 


Ko et al 


as 


uids in the CR-plasma system (e.g., 


Drury & Volk 


1981 


Axford et al. 


1982 


19971). One may a 


McKenzie & Volk 


so consider the sel 


1982 


Ko 


1991a, b; 


-excited waves 


Ko et al. 


1991 


(see e.g., I 

)rury & Fallc 

1986; Zank & McKenzie 

1987 

Zank et al. 

1990; 

Kane e 

al. 

1992 


Zank 

et al 

1993) 

, and magneto-acoust 

ic instabil 

ity (see e.g., 

McKenzie & Webb 

1984 

; Zank 

1989; 

Ko & Jena 

1994; 

Lo & Ko 2007 

Ko & Lo 

2009) 




The influence of CRs on various instabilities has been studied by means of linear 


magn eto-rotational instability (IKhaienabi 


2014), and also galactic dynamo ([Parker 


2012 


1992; 


Kelvin-Helmholtz insta 


Hanasz et al 


2004 


anal\ 

r sis and 

V1HD simulations, for instance, the evo 

ution of Parker instability ( 

3 arker 

1966 

Hanasz 

1997; 

Hanasz & Lesch 

2000; 

Rvu et al. 

2003 


vuwabara et al. 

2004; 

Lo et al. 

T-1 

t-H 

o 

CM 

. Parker-Jeans instability ( 

Kuznetsc 

v & Ptuskir 

l 1983 

Kuwabara & Ko 

20C 

6), 


jility (S uzu ki et ah 


2009) . The results of 


these works showed that in some cases the growth rate has some intriguing dependence on 
the cosmic ray pressure and the coupling of CR and thermal plasma (i.e., the cosmic ray 
diffusion coefficient). For example, while the cosmic ray pressure may effectively enhance the 
Parker and Kelvin-Helmholtz instabilities, small diffusion coefficient can impede the growth 
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(ISuzuki et al.l l2014r) . Moreover, the diffusion coefficient may determine the fragmentation 
direction of Parker-Jeans instability (Kuwabara & Ko 2006). 


Magneto-rotational instability (MRI) is an impo rtant mechanism in differentially 


rotating astrophy s ical o bjects with magnetic fields. 


Balbu s fe Hawley ( 1991 1 and 


H awle y fe Balbu s (11991 1 showed that local and extremely powerful instability in a 


differentially rotating systems with a weak magnetic held destabilize the systems strongly. 
As MRI occurs in accretion disk, the magnetic energy is amplified inside the disk, and 
angular momentum transfer takes place, which is important for obtaining high enough 
accretion rate to explain observations. The efficiency of angular momentum transport 
can be estimated from the saturation level of the magnetic energy, and 


Sano et ah 


(119981 1 


showed that the saturation level of MRI using the resistive MHD simulations. 


Khajenabi (2012) studied the influence of CRs on MRI in the case of dominant toroidal 


magnetic held in the linear regime, and showed that the CR pressure enhanced the growth 
of MRI and the diffusion of CRs suppressed the growth of MRI. In this work, we analyze 
the case of dominant poloidal magnetic held by linear perturbation analysis and MHD 
numerical simulations. We arrive at a somewhat different conclusion. We find similar 
enhancement of MRI by CR pressure as in Khajenabi ( 2 OI 2 L However, we notice that 


diffusion of CRs enhances the growth of MRI as well. This may be alluded to the fact 
that we are using non-uniform initial equilibrium state and different orientation of the 

in Parker or Parker-Jeans instabilities with CRs 


magnetic held. Similar result is observec 

in P 

(Kuwabara et ah 

O 

O 

CM 

Kuwabara & Ko 

O 

O 

CM 


This paper is organized as follows. In section [2] we describe the two-huid model of 
CR-plasma system. In this section we present the governing equations of the shearing box 
model and the rotating cylinder model, and their equilibrium models. In section [3l the 
linear stability analysis and its results of the two models are presented, and in section [4] the 
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results of MHD simulations are presented, section [5] provides a summary and discussion. 

2. MODELS 

We study the MRI in differentially rotating disk in the context of the two-fluid 
CR-plasma system. CR is considered as a massless fluid but with significant pressure. The 
CR fluid is couple to the other fluid, thermal plasma, through the embedded magnetic 
irregularities or hydromagnetic waves. To a first approximation, the effect of waves is 
contained in the hydrodynamical diffusion coefficient of CR and this diffusion coefficient 
serves as the coupling between the two fluids. The system is governed by the total mass, 
momentum and energy equations for the thermal plasma, cosmic ray and magnetic field. 

The cosmic ray energy equation describes the energy transfer between the plasma and 
CR. In this work, we ignore the cross-field-line diffusion of CRs, as in many cases the ratio 


Giacalone & Jokroii 

1999; 

Rvu et al. 

CO 

o 

o 

CM 


20031) . Moreover, ideal MHD is assumed in this work. 


The cases for cross-field-line diffusion and non-ideal MHD will be considered in subsequent 
work. 


The set of governing equations in rotating frame is: 

| + V.(pV) = 0, 


( 1 ) 


d_ 

dt 


(pV) + V 


pVV + I P g + P c + I BB 


2p 0 


A*o 


9P* 

dt 


dPr 


+ V ■ VP g + 7 g -P g V ■ V = 0 , 


BB 


+ p [2ft xV + flx(Oxr)-g] = 0, 

( 2 ) 
( 3 ) 


+ V • VP C + 7 C P C V ■ V - V • ( k || —— • VP C ) = 0 , 


( 4 ) 


dt 


OB 

dt 


V • (V x B) = 0 , 


( 5 ) 
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where p and V are the plasma density and velocity, P g and P c are the thermal pressure 
and the CR pressure, y g and y c are the polytropic index for the plasma and the CRs (i.e., 
the energy densities of the thermal plasma and CR are given by E t h = P g /(y g — 1) and 
E c = P c /(7c — 1)), B is the magnetic field, B is the magnitude of magnetic held strength, 
n || is the CR diffusion coefficient along the magnetic held, / is the unit tensor, and g is 
gravity and 12 is the angular velocity of the rotating frame. Equation ([5]) is the Faraday’s 
induction equation. The inner product of this equation with B gives the energy equation 
for the magnetic held. 

In the following we adopt two models for the differentially rotating disk: the shearing 
box model and the differentially rotating cylinder model. 


2.1. Shearing box 


We consider a two-dimension shearing box in a rotating frame. We choose the local 
Cartesian coordinates ( x,y,z ), where e x is the radial direction, and the angular velocity of 
the rotating frame is 12 = 12e z (see Figured]). The centrifugal force term together with the 
gravity term in Equation ([2]) is replaced by —2ql} 2 xe x (i.e., put 12 x_[!2 xr) 

This is the tidal expansion of the effective potential (see e.g., 


12 x (12 xr) — 

g = - 

Hawlev et ah 

1995) 


2.1.1. Initial equilibrium state of shearing box model 

We adopt the following state as the initial equilibrium state of the shearing box model. 
Density, plasma pressure, CR pressure, magnetic held strength are taken as constant. The 
components of the magnetic held and the velocity are chosen as 


B x — By — V x — V z — 0 , 


( 6 ) 
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V y = —qltx , 
] 


(7) 



( 8 ) 


where (3 is the initial ratio of the magnetic pressure to the thermal plasma pressure. We 


set the initial CR pressure as P c = aP g and P g = PoC'so/7g * s ^ ie thermal plasma pressure 
(and C s o is the sound speed). Setting q — 3/2 in Equation (J7J) gives the Keplerian rotation. 
We set the units of calculation as follows: the units of density, velocity and length are 
p 0 = 1-6 x 10 -24 g cm -3 , C s0 = 10 6 cm s~ 4 and H 0 = 3 x 10 20 cm, respectively. Figure [3] 
shows the distribution of the normalized physical values stated above. In this example, we 
take 7 g = 5/3, y c = 4/3 and a — 1. 


2.2. Differentially rotating cylinder 


Another model of interest is the differentially rotating cylinder model. In this model 
we use the cylindrical coordinate ( r,(j),z ) and consider the system in inertial frame, i.e., put 
0 = 0 in Equation (J2]) but keeping the gravity term (see Figure |2J). 


2.2.1. Initial equilibrium state of differentially rotating cylinder model 


We adopt the following state as the initial equilibrium in the case of differentially 
rotating cylinder model. The equilibrium distribution of a rotating cylinder is obtained 


from the Newtonian analogue of the relativistic tori of Abramowicz et al.J (1197811 . Since we 


are interested in regions close to the equatorial plane, i.e., z r, hence for simplicity we 
assume that the initial equilibrium state depends on r only, and 


V — + V z e z , B — + B z e z . 


(9) 






We note that in this case the diffusion term in Equation (J3J) vanishes. Momentum balance 


m e r gives 


P g + Pr 

p dr 


(-Bl + B iy 

2p 0 


r V A*o P 


~9r = 0 


( 10 ) 


To illustrate ideas, we take the initial total pressure (sum of thermal pressure and CR 
pressure) in the rotating torus as, 


P g + P c = Psum = Kp 1+1/n , Pc = aP g . 


( 11 ) 


Note that a change in a does not change the density distribution. This is more convenient 
when we analyze the dependence of the MRI growth rate on a. We assume B z = B 0 is 
constant, B^ — 0, g = —VT, and the distribution of specific angular momentum (.L = rV^) 
as 


L = Lr 


r o 


( 12 ) 


then the density distribution of the rotating plasma torus is determined by 

(1 + a)P g 


(n + 1 )- 


Ll 


^1 


+ T =£, 


(13) 


p 2(a — 1 )tq V r o7 

where £ is a constant (cf. Bernoulli theorem in fluid physics). In the rest of the paper, 
we consider the gravitational potential is dominated by a point mass at the center, 

T = — GM/y/r 2 + z 2 ~ —GM/r. We consider a non-rotating high-temperature halo outside 
the rotating plasma torus. We take isothermal equation of state for the halo, and adopt the 
distribution 


P = Ph exp 


'1 


_Ch 

\ r / 


(14) 


where ph is the density of the halo at r = r 0 . Here eh = h/V^cn where C s h and Vko are the 
isothermal sound speed (in the halo) and the Keplerian velocity at r = r 0 . We take r 0 as 
the radius at which the density of the rotating plasma torus is maximum, and this density is 
denoted as po- We set the units of length, velocity, time and density as r 0 , Vko, R)/Vko and 
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Po, respectively. Subsequently, we have two nondimensional parameters for the initial torus 

(15) 


£th — 


L 's0 


Cb = 


V 2 

k ao 

V 2 ’ 

V K0 


7g^K0 ’ 

where C s0 = (bg-Pgo/Po) 1 ^ 2 is the sound speed in the torus at r = r 0 , and Vao = (-Bq/PoPo) 1 / 2 
the Alfven speed at r = tq. y g is the adiabatic index of the thermal plasma in the torus. 
In fact, if we represent the gravitational energy by pV£ 0 /2, then £b is the ratio of magnetic 
energy to gravitational energy at r = r 0 , and e t h is (y g — l)/2 times the ratio of thermal 
energy to gravitational energy at r = r 0 . 

As an example, we pick n — 3 , a — 0 (i.e., L is constant), ph/po = 10 -3 , eh = 1.0, 
e t h = 5.0 x 10~ 2 , 6 b = 4.0 x 10~ 4 . Figure |4] shows the distribution of the normalized 
physical values stated above. In this example, we take y g = 5/3,= 4/3 an d a = 1. The 


equilibrium model presented here is a modification of the one in 
include CRs. 


Kuwabara et ah 


(120051 ) to 


3. LINEAR STABILITY ANALYSIS 

We perform standard linear stability analysis on the set of equations ([I})-©. Recall 
that in the shearing box model the term Q x (f2 x r) — g = —2 qVt 2 xe x , while in the 
differentially rotating cylinder model ft, = 0 and g = — V'f. 

In the following analysis, the unperturbed background we consider depends only on one 
coordinate and the velocity and magnetic field is orthogonal to this coordinate axis (this is 
slightly more general than the initial equilibrium state described in previous section). 
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3.1. Shearing box 

In the shearing box model, we denote the set of physical quantities of 
interest as x = {p,V x ,V y ,V z , P g , P c , B x , B y , B z } and the perturbed quantities 
S\ = {Sp, SV X , SV y , 6V Z , 6P g , 6P C , 6B X , 5B y , 5B Z }. We consider the perturbation of the 
form 

Sx(t, x, y, z ) = 5x(x) exp (at + i k y y + i k z z ) , 


where dx = {dp, dV x , i dV y , i dV z , dP g , dP c , —i 5B X , SB y , 5B Z }. After some manipulations, the 
set of linear perturbation equations can be reduced to two first order ODEs. In fact, these 
two ODEs are the continuity equation and the momentum equation. Explicitly, 


where 


d SV X 


An A 1 
Aoi An 


SPt = dP g + dP c H- (.BySBy + B z dB z ) , 

ho 


. AW? B v ( kyB v + k z B z 
1 (1 + W)E Jb* + B 2 ) 


(1 + W)A 2 [ p dx (1 + W )E (B 
2 Vtk v i d 

+! (1 + WOE + E Tx (W + KVz) ' 


A 12 — 


g _ {k 2 y + g) 

p(l + W) 2 A 2 p(l + W)E 


(l + W)E 2 + 20- 


4Q 2 1 rfpdP t 

(1 + W) p 2 dx dx 


. 2klV'jB y (k y B y + k z B z ] 
* (1 + IE)E(P2 + P2) 


A22 — 


1 1 dP t . %W}B y (k y B y + k z B z ) . 20/q, 

(1 + W)A 2 ~~p~d^~ % (l + W)Z(B 2 + B 2 ) ~ 1 (1 + W )E ’ 
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£ = cr + i kyVy + i k z V z , 

(23) 

n 2 T/2 

A 2 - C 2 + c + A 

s (1 + K) (1 + W)’ 

(24) 

VI (k,B, + k,B,f 
(Bl + BJ) ’ 

(25) 

^|| ^zBz) 

E (Bl + Bl) ’ 

(26) 

7g-Pg Z~l2 7c Pc ^ r 2 {By + B z ^) 

P P Pop 

(27) 

{Bl + Bl) 

P t = P g + P c + - z -l . 

ZflQ 

(28) 


The other perturbed quantities can be expressed algebraically in terms of 5V X and 5P t (see 
Appendix A). 


3.1.1. Result of shearing box model 

We take the initial equilibrium state described in section 12.1.11 as the unperturbed 
state, the boundary conditions at x — 0.25 in Figure [3] are taken as 5V X = 1 + iO and 
5P t = 0 + i 0. This condition allows perturbation of the flow to pass through the boundary 
in the x-direction. Moreover, the total pressure is held constant on this boundary. On the 
other boundary at x = —0.25, we require 9ft(<514) 0 and 3ft(AP t ) = 0. (This carries the 

same meaning as the conditions at x = 0.25.) 

We solve the set of linearized perturbation equations, the set of ODEs by shooting 
method. For a trial value of a, we integrate each equation from the boundary at x = 0.25 
(with the assigned boundary value) to the boundary at x = —0.25. We then adjust the 
value of a until 5P t matches the boundary condition at x = —0.25. We take this value of 
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a as the eigenvalue, and take the maximum value of a as the maximum growth rate of the 
system. 


Figure [5] shows the result of the linear stability analysis of the shearing box model. 
The figure displays the dispersion relation for different CR diffusion coefficient km. In the 
figure, cr is the growth rate, and k z is the wave number in the the direction of the initial 
magnetic held. Here we take the CR diffusion coefficient as an input parameter and other 
quantities as fixed parameters (e.g., the ratio of the CR pressure to the gas pressure a — 1, 
the ratio of the gas pressure to the magnetic pressure /3 = 100, the rotational angular 


frequency Q — 1). The maximum value of the normalized n 


to K | 


2001 


= 3 x 10 28 cm 2 s , the value estimated in our Galaxy (IBerezinskii et al. 


= 200 in Figure d corresponds 


1990 


Ptuskin 


Ryu et ah 


2 003 ). The maximum growth rate is given at k z ~ 8.8 and the cut-off 


wave number where the growth rate becomes zero is k z ~ 15.9. In Figure |5l The dispersion 
relations for different k\\ almost completely overlap each other, therefore, we can see only 
one curve in this scale. Figure [6] shows the dispersion relation for different a. In this figure, 
the value of km = 200 is fixed and the other parameters are the same as in Figure El The 
dispersion relations for different a also almost completely overlap each other. We point out 
that the two profiles of Figures El & El are the same. In this model, neither the ratio of CR 
pressure to thermal pressure (while the sum is kept constant) nor the diffusion of CR will 
affect the growth rate significantly. 


3.2. Differentially rotating cylinder 

In the differentially rotating cylinder model, we denote the set of physical quantities 
of interest as x' — { P , V r , V^, V z , P g , P c , B r , B z } and the perturbed quantities 
5x' = {5p, 6V r , 8V$, 5V z , SP g , 5P C1 SB r , SB^, 5B Z }. We consider the perturbation of the form 

Sx'(t, x, y, z ) = Sx'(x ) exp (at + i m<j) + i k z z ) , (29) 
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where 8x' = {8p,8V r ,i8V 4> ,i8V z ,8P g ,8P c ,—i8B r ,8B 4) ,8B z }. Again the set of linear 
perturbation equations can be reduced to two first order ODEs, and these two ODEs are 
the continuity equation and the r-momentum equation. Explicitly, 


d 

8V r 



A' 

A V2 


dK 

dr 

1 

< 

°o 


^21 

A' 

^22 


-1 

< 

°o 


where 

8Pt = 8P g + 8P C + — (B^Bf + B Z 8B Z ) , 
ho 


(30) 


(31) 


A' n = 


(1 + W')A' 2 


i dp; , (i - w')vp b\ 


+ 


2 VtVl 2 B, 


—i 


p dr (1 + W')r (. B 2 + B 2 ) 

A 


(1 + W')iy (B 2 + B 2 ) V r 
2 V^mB^ rm 


( — B& + k z B z 


+i 


(1 + hE')E ,2 r 2 (B 2 + B 2 ) V r 
2mQ 1 | 1 dE' 
(1 + W’)Y,'r ~ r + T/~dr ’ 


-B..;, + 


(32) 


^ _+ 

12 p(l + W') 2 A' 2 p{ 1 + fE')S' V r 2 2 ' ’ 


(33) 


4' 

^21 — 


— ((1 + fE')E' 2 + 2rQ— + —- 

E' ; dr (1 + W') 


4 

r 


(1 - W')VpBl 
(1 + W')r (B: - B)| ' (1 + H’')E' (BJ + B)J V r 


— I 


2W 2 b 4 


7 TI 

—+ k z B : 


v{ 2 b 2 


r {B 2 + A?) 


2 dB# ldp 2(1-W') 


1 

~J\f 2 


1 dPi 


+ 


B ^ dr p dr r(l + hh') 

(1 - W")K 2 B| 


+ 


1 dp dP( 
p 2 dr dr 


p dr (1 + W’)r (B 2 + B 2 ) 


—i 


2 nv^ 2 B 4 


(1 + W) E' (B 2 + 5 2 ) V r 


/ 777 , 

( — Pa + 13-Bj 


(34) 
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A' 

^22 — 


(1 + W')A' 2 


1 dP! 


(1 - W')V' 2 B\ 


+ 


P dr (1 + W')r [B\ + B*) 


2 nv{ 2 B 4 


m 


(1 + W')V (B* + Bl) V r 

2 V'p^mB^ 


— B 0 + k z B : 


(1 + W')YJ 2 r 2 ( B | + Bl) V r 


/ yyi 

( — B (h + AxIT ) + i 


2mBl 


(1 + W')H'r 


(35) 


and 


The other 


S' = cr + * mfl + i k z V z , 

(36) 

a*-c* + c < + A 

s (1 + IP) (1 + W) ’ 

(37) 

W'= . ( Bj, + k z B z ) , 

E' 2 (Bl + Bl) V r 0 V ’ 

(38) 

K = —t—--- ( —Bs + ) , 

E' (Bl + Bl) \r 0 V ’ 

(39) 

Tg^g 7oft „,2 (BJ + B 2 ) 

— 5 — ? V a — ? 

P P Pop 

(40) 

p, P , r , W + B-) 

1 g + c + 2/i 0 ’ 

(41) 

n = 5. 

(42) 


r 


perturbed quantities can be expressed algebraically in terms of 5V r and 5P t ' (see 


Appendix A). 


3.2.1. Result of Differentially rotating cylinder model 

We take the initial equilibrium state described in section 12.2.11 as the unperturbed 
state. Similar to the shearing box model, the boundary conditions at the outer boundary 
r = 4.0 in Figure [I] are taken as 5V r — 1 + % 0 and 8P t ' — 0 + i 0. Hence perturbation of the 
flow can pass through the boundary in the r-direction. Moreover, the total pressure is held 
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constant on the outer boundary. At the inner boundary r = 0.4, we also require 9?(5W) % 0 
and 3f 1(5P t ') = 0. 


Similar to the case of the shearing box model, we solve the set of linearized perturbation 
equations of the differentially rotating cylinder model Equation 1301) by shooting method. 
For a trial value of cr, we integrate each equation from the boundary at r = 4.0 (with the 
assigned boundary value) to the boundary at r = 0.4. We then adjust the value of a until 
SP t ' matches the boundary condition at r = 0.4. We take this value of o as the eigenvalue, 
and take the maximum value of a as the maximum growth rate of the system. 


Figure [7] shows the result of the linear stability analysis of the differentially rotating 
cylinder model. The left panel of the figure displays the dispersion relation for different 
CR diffusion coefficient k\\. Here /cy = 0.4 corresponds to the nominal value in our Galaxy 
K\\ = 3 x 10 28 cm 2 s -1 . The cut-off wave number where the growth rate becomes zero takes 
the same value for different values of /cy except when /cy = 0.0. I 11 the case of k» = 0.0, the 
cut-off wave number is about 5.4% smaller. This can be traced back to the fact that the 
unstable mode of the non-zero /cy case (between the two cut-off wavenumbers) becomes 
neutrally stable (growth rate equals zero) when /cm turns to zero exactly. The cut-off wave 
number in the case of /C|| = 0.0 is smaller because the unstable criterion depends on the 
combine pressures of plasma and CRs (compare to plasma pressure only in th e case of 


/c y > 0.0 as CR diffuse through the plasma). Similar result was obtained in 


Kuwabara fe Ko 


(120061 ) for the role of CRs on Parker-Jeans instability. The maximum growth rate becomes 
larger as /cy increases. The right panel of Figure [7] shows the dependence of the maximum 
growth rate on /C||. Note the horizontal axis is in log scale. The maximum growth rate 
does not change much when /cy < 0.0005, then it increases considerably in the range 
0.0005 < /C|| < 0.05, and then kind of saturated when /cy > 0.05. 


Figure |8] shows the growth rate dependence on a, the ratio of CR pressure to thermal 
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plasma pressure. In this figure, the diffusion coefficient is fixed at ku = 200. The larger is a 
the larger is the growth rate and the larger the cut-off wavenumber. 

4. 2.5-DIMENSIONAL SIMULATION 


In this section, we solve the MHD equations combined with the CR energy equation, 
Equations (JT|)— (j5j) , by MHD simulation code augmented with CR. For the shearing box 
model we put the term fi x (ff x r) - g = —2 qQ 2 xe xl and for the differentially rotating 
cylinder model we set S2 = 0 and g = —VT. The MHD simulations are 2.5-dimensional 
nonlinear, time-dependent, and compressible in cartesian coordinate for the shearing box 


model, and in cy 


Kuwabara et al 


indrical coordinate for the differentially rotating cylinder model. In 


(2004), we used a hybrid scheme to simulate the CR-MHD system. We 


used the Lax-Wendroff scheme for the MHD part and the biconjugate gradients stabilized 
(BiCGstab) method for the diffusion part of the CR energy equation as described in 


Yokovama fe Shibatal (120011 1 to reduce computation time. However, in this work we use the 


Lax-Wendroff scheme for all the equations (MHD and CR equations), because computer is 
very powerful nowadays. The calculation time for such 2.5-dimensional simulation is rather 
short. 


We adopt the MR 

3 code developec 

by 

Hiibata 

(1983 

) and subsequently extended by 

Matsumoto et ah 

(1996 

); 

Havashi et al. 

(199( 

3). Currently, this MHD code is incorporated 


in the Coordinated Astronomical Numerical Software (CANSjj and anyone can use it under 
the acceptance of their licenses. 


Tittp://www. astro.phys.s.chiba-u.ac.jp/cans 
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4.1. Numerical results of the shearing box model 

In the shearing box model, we calculate within the region extracted from x-z plane as 
shown in Figure [D The size of this region is 0.5H 0 x 1.0H Ol with x € [—0.25.f/o, 0.25// 0 ] 
and z G [O.Of/ 0 , 1 -0/T o ]. The numerical grid resolution and the grid size are N x = 41, 

N z = 82, and Ax = 0.0125// 0 , Az = 0.0125f/o- We assume a periodic boundary at 
x = x m ; n , x = x max , and at z — z m j n , z = z max . The initial equilibrium state is described in 
section 12.1.11 To start the simulation, a small velocity perturbation is added to the initial 
equilibrium as follows, 

SV y = 10 -3 x nin(k z z) . (43) 

We choose k z = 10 as a reference to the result of linear analysis (see Figure EJ). 

We study two values of the CR diffusion coefficient, k\\ = 10~ 4 , and k\\ = 10.0 as the 
representative values (see the right panel of Figure [ 3 ). We should point out that Figured 
is the result of linear stability analysis of the differentially rotating cylinder model. The 
maximum growth rate <r max is low for k\\ = 10 -4 , while it is high for /cy = 10.0. In fact, the 
linear analysis on the shearing box model showed that the growth rate is almost the same 
for different /cm (see Figure [5]). This is confirmed by MHD simulations (see below). 

Figure [9] shows the time evolution of the distributions of the magnetic field and the 
CR pressure. In the figure the white curves are the magnetic field lines and the gray-scale 
contour shows the CR pressure. The top three panels show the time evolution for the case 
of K|| = 10 -4 , and the bottom three panels for the case of k\\ = 10.0. The time evolution of 
the magnetic field lines looks like almost the same even if the values of Ku are different. On 
the other hand, the CR pressure distribution are different with different /cm value. In the 
case of /C|| = 10~ 4 , the CR pressure becomes stronger slightly at the valley of the magnetic 
field lines as the time proceeds. However, in the case of ku = 10.0, it shows no variation as 
the time proceeds. 




18 


To compare the results obtained from linear analysis and MHD simulations, we examine 
the temporal variation of V x at a particular point. Figure [TO] shows the time evolution 
of the absolute value \V X \ at (x,z) = (0.0, 0.5). The solid-line corresponds to the case of 
/C|l = 10 -4 , the dash-line corresponds to the case of fty = 10.0, the dotted-line correspond to 
the power-law relation given by the linear analysis. The solid-line and the dash-line almost 
completely overlap with each other and the two lines appear to be one line in this scale. 
The slope of these lines agrees well with the dotted-line from linear analysis. 


4.2. Numerical results of the differentially rotating cylinder model 

In the differentially rotating cylinder model, we calculate within the region extracted 
from r-z plane as shown in Figure [2] The size of this region is I.OHq x 0.5Hq, with 
r G [0.5H 0 ,1.5H 0 ] and z € [O.OH 0 ,0.5H 0 ], The numerical grid resolution and the grid 
size are N r = 81, N z = 42, and Ar = 0.125H 0 , A z = 0.125H 0 . We assume a symmetric 
boundary condition at r = r min , a free boundary condition at r = r max , and a periodic 
boundary condition at z — z min , z = z max . The initial equilibrium state is described in 
section 12.2.11 To start the simulation, a small velocity perturbation is added to the region 
where the rotation velocity is not zero, 

dV^ = — 1CT 3 x cos (k z z ). (44) 

We choose k z = 25.0 as a reference to the result of linear analysis (see left panel of Figure [7]). 
With this choice the analysis of the results of the MHD simulation is easier, because we 
need to control just two waves inside the simulation box. 

Similar to the shearing box model, we also study the two values of the CR diffusion 
coefficient, ku = 10~ 4 , and Kn = 10.0 as the representative values in accordance with the 
result of linear analysis (see the right panel of Figure [7]). Figured!] shows the time evolution 
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of the distributions of the magnetic held and the CR pressure. In the figure the white 
curves are the magnetic held lines and the gray-scale contour shows the CR pressure. The 
top three panels show the time evolution for the case of Ku = 10~ 4 , and the bottom three 
panels for the case of /?u = 10.0. 

In the case of small diffusion coefficient ku = 10~ 4 , the growth of the instability is 
slow. It is still rather insignificant around t ~ 2.0, and the instability starts to grow around 
t = 3.0 (see upper panels of Figurc fTOj) . On the other hand, in the case of larger diffusion 
coefficient, the instability is already approaching its saturation around t ~ 3.0 (lower 
panels of Figure fTTl) . As the growth of the instability proceeds, the low CR-pressure region 
penetrates into the high CR-pressure region around t ~ 3.0. 

In order to understand the mechanism causing different growth rate of MRI, we 
compared the case of ku = 10.0 with the case of /cm = 0.01. They show similar growth 
process of the instability in magnetic fields except that the growth rates are different. The 
left panels of Figure fl2l shows the density (gray scale contour), velocity distribution (white 
arrows), and a reference magnetic held line (white curve) for ku = 0.01 at t — 3.45 and 
/C|| = 10.0 at t — 3.0. The black arrow at the top-right corner is half the unit velocity, 
the Keplerian rotation speed at r = 1.0. High density region is created where the MRI is 
growing strongly. The right panels of Figure [T2] shows the CR pressure distribution, the 
density distribution, and the toroidal velocity distribution along a reference magnetic held 
line for k\\ = 0.01 and 10.0. The CR pressure distribution differs significantly for different 
Ku. For large Ku the CR pressure becomes uniform along the magnetic held line, while 
for small Ku the CR pressure varies in sync with the plasma density. Density attains its 
maximum at the region where the MRI is growing strongly, and its value is higher for the 
larger k\\. The toroidal velocity varies anti-sync with density, but the distributions for 
different k\\ are more or less the same. 
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5. SUMMARY AND DISCUSSION 


We studied the MRI with the effect of CRs by linear stability analysis and MHD 
simulation. We examined two different models: the shearing box model and the differentially 
rotating cylinder model. 


In linear stability analysis, we reduced the set of perturbation equation to two first 
order ODEs and obtain the dispersion relation using shooting method. For the shearing 
box model, the growth rate barely depends on the value of /t y (see Figure [5]). This i s star kly 


different from previous stud ies on related topics (e.g., 


2004 


Kuwabara fe Ko 


Ryu et al 


2003 


Kuwabara et ah 


20061) . which showed considerably dependence of the growth rate 


on the value of /cy. The reason lies in the distribution of CR pressure distribution in 
the initial unperturbed background. If the CR pressure is uniform distributed in the 
unperturbed background (as in the case of the shearing box model), then the growth rate 
will be (almost) independent of the value of Ky. However, for non-uniform CR pressure 
distribution, the growth rate will depends on ku. We confirmed this in our second model, the 
differentially rotating cylinder model, which has a non-uniform CR pressure distribution in 
the unperturbed background. Figure [7] shows the dependence of the growth rate on /cy. The 
growth rate increases as /cy increases, and saturated at large k\\ (see right panel of Figured 
for the maximum growth rat e). This is consistent with the studies on Pa rker instability 


and Parker-Jeans instability (Kuwabara et al. 


2004 


Kuwabara fe Ko 


2006]). However, there 


are some subtle differences. At small values of /cy (< 0.001), the maximum growth rate is 
more or less the same in MRI (see right panel of Figured), but this characteristics was 


not observed in the study of Parker instability ( Kuwabara et al. 


2004 ). Figure |8] shows 


the dependence of the growth rate on the ratio of CR pressure to thermal pressure a. 
The growth rate increases as a increases samely in «y. An increase of a is equivalent to 
a decrease of the ratio of thermal pressure to magnetic pressure. This result is somewhat 
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that the growth rate becomes larger as the 
ratio of thermal pressure to magnetic pressure is larger. This difference is perhaps come 
from our formalism and the non-uniformity of the unperturbed state. In our treatment (see 
Equation [II]) the density distribution is independent of a once we keep the the sum of CR 
pressure and thermal pressure fixed. It is more convenient to study the effect of a without 
changing the density profile. 

In the MHD simulation for the shearing box model, we also obtained the result that 
the growth rate of MRI does not depend on the /cm (see Figure E]). In Figure [TUI we 
compared the growth rate obtained from the linear analysis with that obtained from the 
MHD simulation, and they agreed well. From these results (linear analysis and MHD 
simulations), we can conclude that the growth of the MRI does not depend on the value 
of the CR diffusion coefficient /cy when the initial background CR pressure distribution is 
uniform, at least in the linearly growing phase. 

In the MHD simulation for the differentially rotating cylinder model, we fold that 
the growth rate of MRI under the non-uniform CR pressure background does depend on 
the value of the CR diffusion coefficient /cy. The growth of MRI becomes faster as the 
K\\ becomes larger (see Figure [TIT). This result is consistent with that obtained from the 
linear stability analysis. This result shows that the MRI with cosmic-ray diffusion strongly 
depends on the distribution of the CR pressure background. If the distribution of CRs is 
non-uniform, the growth rate of MRI may change drastically with the value of /cm. 

In the differentially rotating cylinder model, the dependence of the MRI growth rate 
on the value of /cy is caused by the difference in CR pressure distribution along a magnetic 
held line. A general property of diffusion is to smooth out irregularities and to reduce the 
gradient of the relevant quantity. If the diffusion coefficient is large (i.e., weak coupling 
between plasma and CR), the CR pressure (or CR energy) approaches uniform distribution 


different from the result by ( Khaienab 
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22 


quickly even if it were driven away from uniformity by the growth of MRI. Under such 
circumstances, the CR pressure gradient along a magnetic field line becomes small and is 
not able to curb the outward movement of plasma by the centrifugal force. Consequently, 
high density region is formed at the location where MRI is growing and the magnetic field 
line develops the loop like structure. If the diffusion coefficient is small, the CR pressure 
maintains non-uniformity longer and hinders the outward movement of the plasma. Hence 
the density is smaller at the location where MRI is growing when compare with the large 
diffusion coefficient case. On the other hand, the toroidal velocity distribution is not 
sensitive to the value of k» (see right panels of Figure [T2j) . This means that the depicted 
magnetic field line in the case of small or large diffusion coefficient (/C||=0.01 or 10.0) rotates 
with the same rotation speed profile. Therefore, the centrifugal force becomes larger at the 
higher density region and the growth rate becomes larger. 

From these results, we speculate that the effect of CRs on MRI will be weak in the 
phase that the turbulence is sufficiently grown up and the distribution of CR pressure 
approaches uniform. Only in the phase when the turbulence is still growing and the CR 
pressure is non-uniform will the effect of CRs on MRI become significant. 

CMK is supported, in part, by the Taiwan Ministry of Science and Technology grant 
MOST 102-2112-M-008-019-MY3. 

A. Perturbation quantities 

As mentioned in the main text, the set of perturbation equations can be reduced to 
two first order ODEs of 5V X and SP t in the case of shearing box model, and SV r and SPt in 
the case of differentially rotating cylinder model. The other quantities are related to these 
two quantities algebraically. We list them here explicitly. 
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A.l. Shearing box 

First, we express dP g , 5P C , SB X , SB y , SB Z , SV y and SV Z in terms of 6V X , SP t and dp, 
then Sp in terms of SV X and dP t . 


iA-If C f ^ - 

dP g 

E \ dx 

dx / 

1 r <4 2 dp 

dP c 

E |_(1 + il) dx 

dx 


~i2 dP 


SV x + pC , 
p 


sv x + 


pCl Sp 
(1 + K) P ’ 


(Al) 


(A2) 


SB X = - — ( k y B y + k z B z ) SV X , 


(A3) 


sB ” = (1 + W)E \ 5 k » (k ’ A + kA) — + sa dp 


p 


p 


B 


1 dp (1 + W) dB y 


SB Z = 


y 1 p dx 

1 


B y dx 
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_ 1 -£T (k y B y + k z B z ) } SV X 


f 1 , „ „ . , „ N 5 -ft , 


x^k z ( k y By + k z B z ) -h Si? 


(1 + 1F)E (E 

'1 dp (1 + W) dB 
p dx 


+B Z 


B z dx 


P 

*1 <514 


P 


(A4) 


(AS) 


did, = 


1 / , SP t VlBy(kyBy+k Z B Z )Sp 

(1 + WOE \ y P (By + Bl) p 

VlByjkyBy + k.B,) d^ dVy 

E (52 + 52 ) p dx v ; dx 


<514 


(A6) 


<514 = 


r 


,, dP t , F A 2 P 2 (fcyP* + **P*) dp 

1 / . ^o\ 

P 


_ y _u _ 1 1 A ^ y y 

(1 + WOE \ z p (Bl + Bl) 


VaBz ( kyB y + k z B z ) dp + ^ d!4 


S (Bl + Bl) p dx 


dx 


SVr 


(A7) 



and finally, 


P 


If 1 6 P t 1 ridP t 

A? \(1 + W) p + E [p-fa 
2 QV^By ( k y By + k z B z ) 
+l {l + W)Y,{Bl + Bl) 


A 2 dp 

p dx 
SV X \ . 


The other quantities are given by Equations flTTj) - (l28]) . 


(A8) 


A.2. Differentially rotating cylinder 


Similarly, we express SP gl SP C , 8B r , SB^, SB Z , SV^ and 5V. in terms of 5V r , 5P t and 
5p, then Sp in terms of 5V X and 5P t '. 


1 (^dp dP s \ sVr+ ,Sp 

P 


SP s=»[ C -Tr ir 


(A9) 


SPr = ~ 


C 2 dp dP c 


S' L(l + K') dr dr 


5V r + 


PCI Sp 
(1 + IP) p ’ 


(A10) 


— 1 / 777 

8B r = --[-B^ + k z B z )SV ri 


(All) 


SB, = 


m fm „ , „ \ SP t , Sp 


(1 + W')Ti' I E'r V r 

a 


( — B, + k z B z 


+ E 'B^ 
P P 


l dp (1 + W') dB$ (1 + W') 

p dr Bfj, dr r 
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p dr B z dr 


SV r 
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<514 = 


1 j m6ID t Vj 2 B$ (m , kR \Sp 

(1 + 14')S' | r p (52 + 52 ) V r 0 2 V p 


/m \ f l dp 2 

S' (Pj + 52 ) v7 0+ 2 VVpd^ + r 


+i(l + T4')r^ + i2n 

dr 


<54. 


(A14) 


<514 = 


(1 + W) £' 




P (S| + Bf) V r 


n 2 B. 


wjBffkfTp^ (7 B * +kA ) + ' (1 + ^ 


dK 

dr 


<514 


(A15) 


and finally, 


<5p _ 1 [ 1 <5P t ' 1 

7 ~~ 77(1 + 14 ') P + 7 


ldP t ' ^' 2 dp (1 — 14') 47-00 


p dr p dr (1 + 14') r ( B | + Bfj 


+i 


2nv{ 2 B 4 


(1 + 14') S' (P| + P 2 ) v r 


771 


Pa, + L-B. 


<5K 


(A16) 


The other quantities are given by Equations (130|) - (I4T]) . 
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Normalized value 



Fig. 2.— Initial distribution of the physical quantities in the shearing box model. 
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Distribution of the physical values 



Fig. 4.— Initial distribution of the physical quantities in the differentially rotating cylinder 


model. 















Fig. 5.— Dispersion relation for the magnetorotational instability with the effect of CRs 
for different /cy in the shearing box model. Here a is the growth rate of perturbation and 
k z is the wavenumber along the direction of the magnetic held in the unperturbed state. 
Apparently, all the cases collapse to one line. This indicates the dispersion relation is almost 
independent of /cy in the shearing box model. The reason is CR pressure is uniform in the 
initial unperturbed state. 
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Fig. 6.— Same as Figure [5] except that the CR diffusion coefficient is fixed at /ty = 200 and 
the ratio of CR pressure to thermal pressure a varies (while the sum of them is kept fixed). 
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kz K|| 


Fig. 7.— Left: Dispersion relation for the magnetorotational instability with the effect of 
CRs at different k\\ in the differentially rotating cylinder model. Here a is the growth rate 
of perturbation and k z is the wavenumber along the direction of the magnetic field in the 
unperturbed state. Right: Dependence of the maximum growth rate <7 max on /cy. <r max is 
small when /cy < 0.0005. it increases considerably in the range 0.0005 < /cy < 0.05, and goes 
to saturation when /cy > 0.05. 
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Fig. 8.— Same as the right panel of Figure [7] except that the CR diffusion coefficient is 
fixed at k\\ = 200 and the ratio of CR pressure to thermal pressure a varies (while the sum 
of them is kept fixed). 
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K||=10.0 


Fig. 9.— Time evolution of the CR pressure distribution and magnetic field lines of the 
shearing box model for the cases of k,\\ = 1CT 4 (top) and 10.0 (bottom). The gray scale and 
the white curves show the CR pressure distribution and magnetic held lines, respectively. 
The magnetic held lines behave almost the same in both cases. However, as time proceeds, 
the CR pressure becomes slightly larger at the valley of the magnetic held lines for the small 


K|| case. 
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Fig. 10.— Comparison of the growth rate obtained from MHD simulations and linear 
analysis in the shearing box model. The results of the two simulations almost overlap each 
other completely (it is difficult to distinguish them in this scale). This verifies the conclusion 
of the result of the linear analysis (see Figure [5]). The line exp[0.75 * (t — 10.0)] shows the 
power-law relation given by the linear analysis. The result of MHD simulations and linear 
analysis agree well with each other. 
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Fig. 11.— Time evolution of the CR pressure distribution and magnetic field lines of the 
differentially rotating cylinder model for the cases of k\\ = 10~ 4 (top) and 10.0 (bottom). The 
gray scale and the white curves show the CR pressure distribution and magnetic held lines, 
respectively. The growth of the instability is slower in the case of small diffusion coefficient 
when compare to the large diffusion coefficient one. 





































































































































